Confined Cubic Blue Phases under Shear 
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Abstract. We study the behaviour of confined cubic blue phases under shear flow via 
lattice Boltzmann simulations. We focus on the two experimentally observed phases, 
blue phase I and blue phase II. The disclination network of blue phase II continuously 
breaks and reforms under shear, leading to an oscillatory stress response in time. The 
oscillations are only regular for very thin samples. For thicker samples, the shear 
leads to a "stick-slip" motion of part of the network along the vorticity direction. 
Blue phase I responds very differently: its defect network undergoes seemingly chaotic 
rearrangements under shear, irrespective of system size. 

PACS numbers: 47.11. Qr, 47.57.Lj, 47.57.Qk, 61.30.Dk, 61.30.Pq, 61.30.Jf, 83.80.Xz 
1. Introduction 

Cholesterics are liquid crystals in which the local nematic director field spontaneously 
twists in thermodynamic equilibrium [lj. The preferred configuration close to the 
isotropic boundary features twist around two perpendicular axes 5 as opposed to just 
one axis in the regular cholesteric state, and the corresponding deformation is denoted 
a "double-twist cylinder". As it is topologically impossible to cover continuously 3D 
space with double-twist cylinders, defects arise. These organise into a variety of regular 
periodic lattices, giving rise to the so-called cubic blue phases (BPs) [2j [3]. There are 
two experimentally observed cubic blue phases, BPI and BPII (a third, BPIII, is thought 
to be amorphous [1]). 

BPs have been long considered as purely of academic interest due to their very 
narrow range of stability. This view has changed since the creation of polymer- 
stabilised and temperature-stabilised BPs [5j [6], which has opened up the possibility 
of novel applications. During the last few years considerable progress has been achieved 
regarding the behaviour of BPs in confined geometries [Zl El E], under external fields 
pLOl [TT1 [T2l [T3l [2] , and in the presence of colloidal particles [TBJ . The kinetics of BP 
domain growth have been recently addressed in [16]. However, our understanding of 
their dynamical behaviour under flow remains very limited. The aim of this work is to 
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address this issue by studying, for the first time, the response of confined BP samples 
to a shear flow. 

Flow response in cholesterics is both strongly non-Newtonian and highly 
anisotropic. For example, if a cholesteric helix is subjected to a Poiseuille flow along 
the helical axis, small pressure differences drive flow mainly through "permeation" , first 
investigated by Helfrich [TT] . In the permeation mode the liquid crystal flows while 
leaving the director field virtually unchanged, which leads to high dissipation and large 
viscosities. Marenduzzo et al. [181 HH] simulated shear and Poiseuille flow in cholesteric 
liquid crystals in the permeation mode, and showed the importance of the boundary 
conditions in determining the apparent viscosity of the fluid. They also found that a 
strong secondary flow appears. Rey [20j [21] studied shear in cholesterics oriented with 
the helix along the vorticity axis and found that, at low Ericksen number, travelling 
twist waves appear which lead to the rotation of the cholesteric helix. At higher forcing, 
the helix uncoils and leaves a flow-induced nematic phase. Rey also studied cholesterics 
subjected to both steady flow and low frequency small amplitude oscillatory shear for 
different helix orientations [221 [23]. He found that splay /bend/twist deformations were 
excited when the helix was aligned along the flow direction; splay/bend deformation 
occurred when the helix was aligned along the velocity gradient; but only twist 
deformations appeared when the helix was aligned along the vorticity axis. 

Dupuis et al. [24] performed the first numerical investigation of BP rheology in 
Poiseuille flow, starting from equilibrium structures of BPI and BPII and a periodic 
array of doubly twisted cylinders. Under small forcing, the network opposed the flow 
giving rise to a significant increase in apparent viscosity. Upon increasing the forcing 
they found clear evidence of shear thinning. In the crossover region they predicted a 
novel oscillatory regime where the network continuously breaks and reforms as portions 
of the disclinations in the center of the channel move to neighboring cells and relink 
with the parts of the network left behind by the flow. Compared with the cholesteric 
case, the viscosity still decreases with forcing (the system shear thins) but much less 
than for cholesterics in the permeation mode, which is in agreement with experiments 
[251 [26]. 

The purpose of the present paper is to contribute to the understanding of the 
flow behaviour of cholesteric blue phases by reporting results of the first supra-unit cell 
simulations. This paper is organised as follows. In Sec. [2] we introduce the Beris-Edwards 
model for hydrodynamics of liquid crystals, which is a generalisation of Ericksen's and 
Leslie's theory of nematodynamics pQ. In Sec. [3] we describe the simulation methodology 
and the boundary conditions we applied. Results for a confined BPII in rectilinear shear 
flow are presented in Sec. |3.1| with particular focus on the dependence of the flow on the 
thickness of the sample and thermodynamic state. In Sec. |3.2| we show corresponding 
results for BPI, which indicate fundamental differences in the flow behaviour of BPI and 
BPII. 
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2. Model 

Our approach is based on the well-established Beris-Edwards model for hydrodynamics 
of cholesteric liquid crystals [27J, which describes the ordered state in terms of a traceless, 
symmetric tensor order parameter Q(r). In the uniaxial approximation, the order 
parameter is given by Q a p = q s {n a n^ — | 8 a p) with n the director field and q s the 
amplitude of nematic order. More generally, the largest eigenvalue of Q, < q s < | 
characterizes the local degree of orient ational order. The thermodynamic properties of 
the liquid crystal are determined by a free energy J 7 , whose density / consists of a bulk 
contribution /& and a gradient part f g , as follows, 

X 2^2 
fg = + ^QoQa/3) + ~^(dpQ a p) . (1) 

The first term contains the bulk-free energy constant A and the inverse temperature 
7 which controls the magnitude of order. The second part quantifies the cost of elastic 
distortions, which are proportional to the elastic constant K\ we work for simplicity in 
the one-elastic constant approximation. The wavevector q = 2n/po, where p is the 
cholesteric pitch. The actual periodicity of the BP structure, p, does not need to be 
equal to p Q . The "redshift" r = p/po is adjusted during the simulation by following a 
simple procedure similar to [28] . 

A particular thermodynamic state is specified by two dimensionless quantities: the 
effective temperature r and chirality ft, which are given by 

27(1 - 7 /3) 

T = 

7 



K 



A)7 

The dynamical evolution of the order parameter is given by the equation 

{dt + u a d a ) Q - S(W, Q) = TH. (2) 

The first term on the left hand side of Eqj2]is a material derivative, which describes the 
rate of change of a quantity moving along with the flow. The second term accounts for 
the rate of change due to local velocity gradients W a p = dpu a , 

S(W,Q) = (£A + n)(Q + ^) 

+ (Q + *)(£A - O) - 2^(Q + \)Tv(QW), (3) 

where Tr denotes the tensorial trace, while A = (W + W T )/2 and ft = (W — W T )/2 
are the symmetric and antisymmetric part of the velocity gradient, respectively. £ is 
a constant depending on the molecular details of the liquid crystal. Flow alignment 
occurs if t; cos 29 = (3q s )/(2 + q 8 ) has a real solution, where 9 is the Leslie-angle: we 
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select this case by setting £ = 0.7 in our simulations. H is the molecular field, which is 
a functional derivative of T that preserves the tracelessness of Q: 

— « + HS- 

The rotational diffusion constant T in Eq. [2] is proportional to the inverse of the 
rotational viscosity 71 = 2q 2 s /Y [I]. 

The time evolution of the fluid density and velocity are respectively governed by 
the continuity equation d t p = —d a {pu a ) ) and the following Navier-Stokes equation: 

d t u a + pupdpUa = dpllap + r]d(3{d a U(3 + dpu a + (1 + 3—)d^u^5 a(3 }. (5) 

The final term in Eq. [5] arises from the Chapman-Enskog expansion of the LB equations 
[29] . At low flow rates the fluid can be considered as incompressible, so that the last 
term on the right hand side of Eq. [5] remains small, r] is an isotropic background viscosity 
which is set to 77 = 0.8333 in LBU. The pressure tensor reads explicitly 

n o j = Po6 a p — iH ai (^Q 1 (3 + 2^7/5^ ~~ £ {Qcx-i + 3^7^ 
/ 1 \ 8J~ 

+ 2£ ( Qa(3 + ^$a(3 J 

+ Qa~fH 7 p — H OL1 Q 1 p. (6) 

In the isotropic state Q = and EqJH] is reduced to the scalar pressure which, in a 
system at rest, is constant to a very good approximation. 

The system of coupled partial differential equations Eqs. [2] and [5] is solved by means 
of a hybrid scheme [30j . This uses a lattice Boltzmann algorithm with predictor-corrector 
scheme for the continuity equation and Eq. [5| and a finite difference scheme for the 
equation of motion of the tensor order parameter Eqj2| More details on the algorithm 
can be found in [29j [3]] . 



3. Methodology and Results 

We present results from simulations of confined cubic blue phases in steady shear flow. 
Our system consists of n x x n z = 4 x 4 unit cells in x- and z-direction (flow and vorticity 
direction respectively) and a varying number n y of unit cells in ^/-direction (which is 
the flow gradient direction in our geometry). As is common in BP simulation studies 
[H [16] , we initialised our simulations with analytical solutions that minimize the free 
energy functional Eq{l]in the high-chirality limit, and equilibrated these (for 3000 LB 
timesteps) with periodic boundary conditions prior to starting the shear. During the 
equilibration run the optimal redshift r was recalculated at every timestep. 

Shear flow was then started by applying shear boundary conditions to the LB 
distribution functions as in p?T] . For Q, the equilibrated order parameter was advected 
on the top and bottom boundary, according to the imposed wall velocities. Periodic 
boundary conditions were retained along the flow and vorticity direction. This set-up is 






Time: 50000 Time: 75000 Time: 100000 

Figure 1. Shear flow of BPII at different gap sizes. The images show snapshots of 
a subset of the entire disclination network in a system with walls at top and bottom, 
and the flow direction into/out of the page. The different number of unit cells in the 
velocity gradient direction have been superposed and are shown in different colours: 
n y = 1 (yellow), n y = 2 (green), n y = 3 (red) and n y = 4 (blue). All cases show 
n x x n z = 2x2 unit cells in flow and vorticity directions. The temperature and 
chirality are r = —0.5 and k = 2, respectively, where BPII is the equlibrium phase. 
The strain rate is 7 = 7.8125 x 10 -5 in reciprocal LB time units in all cases. 



similar to the one Dupuis et al. used to simulate Poiseuille flow [21]. The pinning of the 
defect network at the boundaries may be realised in practice by impurities or surface 
defects. During shear flow the unit cell size was kept fixed and a constant redshift equal 
to the value at the end of the equilibration protocol was assumed. 

The total number of unit cells in our simulations of both BPI and BPII was 
n x x Tty x n z = 4 x {1,2,3,4} x 4. A resolution of 32 lattice Boltzmann units (LBU) 
per pitch length (16 per unit cell) was used for BPII, and 64 LBU per pitch length 
(32 per unit cell) used for BPI. This provided enough spatial resolution to track even 
complex rearrangements of the disclination networks in each case. The timestep and 
lattice spacing in LBU corresponds roughly to ~ Ins and ~ lOnm in SI units. The LBU 
of stress is equal to about 10 8 Pa. More details about the conversion from LBU to SI 
units can be found in [H, [16] . 



3.1. Blue Phase II 

We first report the results for BPII. Fig.[T]shows the disclination network in steady shear 
flow for different numbers of unit cells between the walls. (We prefer the perspective 
along the flow direction because the network topology is clearer.) Temperature and 
chirality are r = —0.5 and n = 2 respectively, and for these parameters BPII is the 
equilibrium phase. The velocity of the walls has been adapted so that the same shear 
rate 7 = 7.8125 x 10 -5 is achieved for all gap sizes. Shear rates four times larger and 
smaller were also applied and produced very similar results, apart from the flow velocities 
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Time: 25000 Time: 75000 Time: 125000 



Figure 2. BPII in shear flow for different thermodynamic parameters. The snapshots 
show the evolution of the network with the flow direction into/out of the page. The 
reference state (blue, r = —0.5, n = 2), a high temperature state (red, r = 0.5, n = 2) 
and a low chirality state (green, r = —0.5, n = 1) are shown; BPII is metastable in the 
later two parameter sets. The strain rate is 7 = 7.8125 x 10 -5 LBU in all three cases. 
The pictures show a central cubic section consisting of 8 unit cells cut from a larger 
simulation with n y = 4. 

which scale up accordingly. For all gap sizes we observe an affine transformation in the 
gradient-flow plane that leads to a break-up of the network at strains of about half a unit 
cell size. Shortly after the break-up the network reforms further downstream, resulting 
in a periodic, regular structure that is very close to the equilibrium BPII disclination 
network. Similar oscillations have been previously reported for BPs in Poiseuille flow, 
with a gap of one unit cell size [24] . 

It is interesting to compare the results for different gap sizes. For n y = 1 we 
observe that the BPII-network is slightly deformed with respect to the quiescent state, 
but always reforms at the same position. For n y > 1, the BPII network distorts more 
under the same shear rate, suggesting that the structure within the sample is strongly 
affected by the anchoring at the boundaries. Interestingly, when the sample is thicker 
than one unit cell, a movement of the network in the positive vorticity direction arises 
which is not observed for n y — 1. This network motion does not proceed at a constant 
speed, but rather proceeds in a "stick-slip" fashion: elastic stresses first slowly build 
up, and are then quickly released as the network "slips" along the vorticity direction. 
Changing from a left-handed to a right-handed helix inverts the direction of motion of 
the network. Hence the "stick-slip" motion is related to the chiral nature of the helix. 
This phenomenon may bear some similarity to the travelling wave motion in cholesterics 
sheared along the vorticity direction below a critical Ericksen number [20l I2T] . 

The speed of the network migration in the vorticity direction depends on the 
thermodynamic state of the blue phase. A comparison between the state from Fig. [j] 
with a "high temperature" state at r — 0.5, ft = 2 and a "low chirality state" at 
r = —0.5, ft = 1 is shown in Fig. [2} Note that in both these states BPII is a metastable 
phase. Fig. [2] shows a sequence of initially overlapping disclination networks which 
separate during the course of the simulation. The network moves fastest in the high 
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Figure 3. BPII: ensemble-averaged shear stress (H xy ) versus time in the reference 
state (r = —0.5, ft = 2) for gaps ranging from one to four unit cells. 

temperature state, is slower in the reference state and slower still in the low chirality 
state. 

The total network displacement along the vorticity direction is more than a unit 
cell in 10 5 timesteps for the high temperature state. Typical average and maximum 
values of the scalar order parameter are q s = 0.338, q s , m ax = 0.383 (low chirality state), 
q s = 0.206, q s ,max = 0.320 (high temperature state) and q s = 0.304, q s ,max = 0.408 
(reference state, in Fig. [I]). The average scalar order parameter of the low chirality 
state is about 10% higher than that of the reference state. The relative differerence 
between the mean and the maximum value is about a factor 2.5 smaller. We propose 
that the speed of network migration is different because permeation, which should be 
required for the stick-slip motion, may occur more easily in states with less order (high 
temperature state) and gentler spatial variation (low chirality state). Furthermore, the 
uneven stick-slip nature of the network motion is possibly due to the anchoring at the 
boundaries which prevents network displacements at the wall. 

Fig. [3] shows the spatially averaged shear stress (H xy ) versus time, in the reference 
state and for different gap sizes. For gaps accommodating only one unit cell a 
regular periodic pattern emerges shortly after a transient start-up phase. Every cycle 
is related to a break-up and reformation of the network. The time-average in the 
steady state is about U xy ~ 5.62 x 10 -5 in LBU and yields an apparent viscosity 
Vapp = (l+n^/777) — 1.86 LBU. This is significantly higher than our baseline Newtonian 
viscosity 77 = 0.8333. When the gap contains n y >2 unit cells, the average stress shows 
less regular oscillations, and fluctuates around a mean value which is about 10 — 15% 
smaller, roughly r] app = 1.65 for n y = 4. The fact that the average stress stabilises for 
n y > 2 indicates that bulk flow behaviour is already reached for those thicknesses, but 
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not in the single unit cell case. Different shear rates lead to similar generic features. 
3.2. Blue Phase I 

We now compare and contrast the results obtained for BPII with those pertaining to 
the rheology of BPI. In equilibrium, the topology of the disclination network in BPI is 
very different from that of BPII. The main difference is the presence in BPII of junction 
points where four disclination lines of topological charge -1/2 meet; no such junctions 
occur in BPI. We give representative results for BPI at r = —0.5 and k — 1, where it 
is the equilibrium phase. Fig. [4] shows the disclination network at timestep t = 4 x 10 5 
for different gap sizes, with the viewing direction again chosen along the flow direction. 
For clarity we only show n x x n y x n z — 1 x {1, 2, 3, 4} x 2 unit cells. 




(b) (c) (d) 

Figure 4. Confined BPI in shear flow at r = —0.5, k = 1 and gap sizes of one (a, 
yellow), two (b, green), three (c, red) and four (d, blue) unit cells. The pictures show 
the defect network along the flow direction at timestep t = 4 x 10 5 at the same flow 
rate as above. 

In contrast to BPII, the disclination network of BPI changes significantly under 
a shear flow, and typical snapshots in Fig. [4] bear no relation to the equilibrium BPI 
configuration. Furthermore, there is no sign of the "stick-slip" motion of the defects 
observed in BPII. While no regular cubic structure emerges under flow, it is interesting 
to note that the network develops a periodic modulation along the vorticity direction 
(horizontal in Fig. [4]). This periodicity arises for all gap sizes. Whether this feature is 
a consequence of frustration caused by the anchoring at the boundaries, or is a prelude 
to a slow ordering of the network under flow, we cannot yet tell. 
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Figure 5. BPI: ensemble- averaged shear stress (IL^) versus time in the reference 
state (r = — 0.5, k — 1) for gaps ranging from one to four unit cells. 

The ensemble averages of the shear stress J\ xy for different gap sizes is shown in 
Fig. [5] and gives more evidence of the chaotic nature of the flow of BPI. The magnitude 
of the average shear stress is comparable to that seen in BPII; the fall in amplitude of 
the oscillations with increasing gap size is also similar. 

4. Conclusions 

In conclusion, we have presented the first supra-unit cell simulations of confined cubic 
blue phases in simple shear flow. The rheological behaviour of BPI and BPII was found 
to be strikingly different. 

Our simulations show that the disclination network of BPII continuously breaks and 
reforms while the network undergoes an affine deformation under shear. This process 
leads to an oscillatory shear stress in time, although the oscillations are only regular for 
gap sizes of a single unit cell across the flow gradient direction, where anchoring affects 
the physics of the system greatly, in agreement with recent simulations of quiescent BP 
networks j7J[8l|9]. For larger gap sizes (two and more unit cells) a "stick-slip" motion 
of the network arises along the vorticity direction. The associated network migration 
is not directly linked to mass flow, and is therefore somewhat akin to permeation 
flow in cholest erics. The average speed of this network migration depends on the 
thermodynamic parameters, and appears to be higher for higher temperature and lower 
chirality. 

In contrast to BPII, BPI shows no regular break-up and relinking of the unit 
cells and its disclination network undergoes chaotic rearrangements under shear. This 
behaviour was found for all gap sizes, and may be linked to the strongly different 
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topology in the equilibrium disclination patterns of BPI and BPII. We will address 
the rheology of BPI in more detail in an upcoming publication [32] . 

We hope that our results will stimulate further experimental work on the rheology of 
blue phases, similar to those reported e.g. in [33] but targeted to micron-sized samples, 
which would provide a direct test of our simulations. 

We finally note that the anchoring we have chosen does not lead to a conflict with 
the bulk ordering of blue phases. In liquid crystal cells, it is customary to treat the 
surface so as to favour normal or planar anchoring of the director field at the boundary. 
When used in very thin samples, these boundary conditions lead to the formation of new 
disclination networks whose rheology it would be interesting to investigate [HJ C2] . 
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